﻿#ifndef JEGEO_MATH_HPP
#define JEGEO_MATH_HPP
#include <algorithm>
#include <cmath>
#include <vector>
#include <iostream>

namespace jegeo::math {

#ifdef JEGEO_REAL_TYPE
using jereal =  JEGEO_REAL_TYPE;
#else
using jereal = double;
#endif


/*  Even though these might be more useful as long doubles, POSIX requires
    that they be double-precision literals.                                   */
#define M_E         2.71828182845904523536028747135266250   /* e              */
#define M_LOG2E     1.44269504088896340735992468100189214   /* log2(e)        */
#define M_LOG10E    0.434294481903251827651128918916605082  /* log10(e)       */
#define M_LN2       0.693147180559945309417232121458176568  /* loge(2)        */
#define M_LN10      2.30258509299404568401799145468436421   /* loge(10)       */
#define M_PI        3.14159265358979323846264338327950288   /* pi             */
#define M_PI_2      1.57079632679489661923132169163975144   /* pi/2           */
#define M_PI_4      0.785398163397448309615660845819875721  /* pi/4           */
#define M_1_PI      0.318309886183790671537767526745028724  /* 1/pi           */
#define M_2_PI      0.636619772367581343075535053490057448  /* 2/pi           */
#define M_2_SQRTPI  1.12837916709551257389615890312154517   /* 2/sqrt(pi)     */
#define M_SQRT2     1.41421356237309504880168872420969808   /* sqrt(2)        */
#define M_SQRT1_2   0.707106781186547524400844362104849039  /* 1/sqrt(2)      */


/**
 * @brief PI相关常量定义
 */
constexpr double PI_DIV_D720 = M_PI / 720.0;
constexpr double PI_DIV_D180 = M_PI / 180.0;
constexpr double D180_DIV_PI = 180.0 / M_PI;
constexpr double PI_MUL_2 = M_PI * 2.0;


/** 公差定义 */
static const jereal EPSILON_HIGH = 1.0E-16;
static const jereal EPSILON_MEDIUM = 1.0E-10;
static const jereal EPSILON_LOW = 1.0E-7;
static const jereal EPSILON = EPSILON_MEDIUM;

template<typename T> T Epsilon();
template<> inline double Epsilon<double>() { return static_cast<double>(EPSILON_MEDIUM); }
template<> inline float Epsilon<float>() { return static_cast<float>(EPSILON_LOW); }
template<> inline long double Epsilon<long double>() { return static_cast<long double>(EPSILON_HIGH);  }
template<> inline int32_t Epsilon<int32_t>() { return 0; }
template<> inline int64_t Epsilon<int64_t>() { return 0; }
template<> inline long Epsilon<long>() { return 0; }

/**
 * @brief 在做数据类型转换时,根据输入输出类型, 进行四舍五入
 */
template<typename T, typename FT>
constexpr inline T roundcast(const FT &v) noexcept { return static_cast<T>(v); }

template<>
constexpr inline int32_t roundcast<int32_t, double>(const double &v) noexcept { return static_cast<int32_t>(v + (v >= 0 ? 0.5 : -0.5)); }

template<>
constexpr inline int64_t roundcast<int64_t, double>(const double &v) noexcept { return static_cast<int64_t>(v + (v >= 0 ? 0.5 : -0.5)); }

template<>
constexpr inline long roundcast<long, double>(const double &v) noexcept { return static_cast<long>(v + (v >= 0 ? 0.5 : -0.5)); }

template<>
constexpr inline int32_t roundcast<int32_t, float>(const float &v) noexcept { return static_cast<int32_t>(v + (v >= 0 ? 0.5 : -0.5)); }

template<>
constexpr inline int64_t roundcast<int64_t, float>(const float &v) noexcept { return static_cast<int64_t>(v + (v >= 0 ? 0.5 : -0.5)); }

template<>
constexpr inline long roundcast<long, float>(const float &v) noexcept { return static_cast<long>(v + (v >= 0 ? 0.5 : -0.5)); }


/**
 * @brief 常用数学函数
 */
template<typename T>
constexpr inline T jeabs(T v) noexcept { return v < 0 ? -v : v; }

template<typename T>
constexpr inline T jemin(T a, T b) noexcept { return a < b ? a : b; }

template<typename T>
constexpr inline T jemax(T a, T b) noexcept { return a > b ? a : b; }

constexpr inline jereal jesqr(const jereal &v) noexcept { return v * v; }

constexpr inline jereal jemul(const jereal &a, const jereal &b) noexcept { return a * b; }

constexpr inline jereal jediv(const jereal &a, const jereal &b) noexcept { return a / b; }

template<typename T=jereal>
constexpr inline T jesqrt(const jereal &v) noexcept { return roundcast<T,jereal>(std::sqrt(v)); }

template<typename T>
constexpr inline T jehypot(const T &x, const T &y) noexcept { return roundcast<T,jereal>(std::hypot(static_cast<jereal>(x), static_cast<jereal>(y))); }


/**
 * @brief 带公差的比较函数
 */
template<typename T>
bool fuzzy_equal(const T &a, const T &b, const T &epsilon = Epsilon<T>()) noexcept {
    T diff = a - b;
    return ((-epsilon <= diff) && (diff <= epsilon));
}

template<typename T>
bool fuzzy_notequal(const T &a, const T &b, const T &epsilon = Epsilon<T>()) noexcept {
    T diff = a - b;
    return ((-epsilon > diff) || (diff > epsilon));
}

template<typename T>
bool fuzzy_greaterequal(const T &a, const T &b, const T &epsilon = Epsilon<T>()) noexcept {
    return (a > b) || fuzzy_equal<T>(a, b, epsilon);
}

template<typename T>
bool fuzzy_lessequal(const T &a, const T &b, const T &epsilon = Epsilon<T>()) noexcept {
    return (a < b) || fuzzy_equal<T>(a, b, epsilon);
}

template<typename T>
constexpr inline bool fuzzy_iszero(const T& v, const T& epsilon = Epsilon<T>()) noexcept {
    return fuzzy_equal<T>(v, static_cast<T>(0), epsilon);
}

template<typename T>
constexpr inline bool fuzzy_notzero(const T& v, const T& epsilon = Epsilon<T>()) noexcept {
    return fuzzy_notequal<T>(v, static_cast<T>(0), epsilon);
}


/**
 * @brief 角度转为弧度
 */
constexpr inline jereal deg2rad(const jereal &deg) noexcept { return deg * PI_DIV_D180;}

/**
 * @brief 弧度转为角度
 */
constexpr inline jereal rad2deg(const jereal &rad) noexcept { return rad * D180_DIV_PI; }

/**
 * @brief 计算圆弧的凸度
 */
inline jereal deg2bulge(const jereal &deg) noexcept {
    // int64_t intdeg = roundcast<int32_t, jereal>(deg * 1.0E8);
    // if (intdeg == int64_t(90e8)) return 0.41421356237309504880168872421;
    // if (intdeg == int64_t(-90e8)) return -0.41421356237309504880168872421;
    // if (intdeg == int64_t(180e8)) return 1.0;
    // if (intdeg == int64_t(-180e8)) return -1.0;
    return std::tan(deg * PI_DIV_D720);
}


template<typename T>
constexpr inline T sign(T v) noexcept { return v < 0 ? -1 : (v > 0 ? 1 : 0); }

template<typename T>
constexpr inline T clamp(T v, T minv, T maxv) noexcept { return jemin(jemax(v, minv), maxv); }


/**
 * @brief 计算圆弧的弧度范围限定在0~2PI之间
 */
inline jereal normalize_radian(jereal rad) noexcept
{
    while (rad < 0) {
        rad += 2 * M_PI;
    }
    while (rad >= 2 * M_PI) {
        rad -= 2 * M_PI;
    }
    return rad;
}

inline jereal normalize_degree(jereal deg) noexcept
{
    while (deg < 0) {
        deg += 360.0;
    }
    while (deg >= 360.0) {
        deg -= 360.0;
    }
    return deg;
}

/**
 * @brief 判断弧度是否在两个弧度之间
 */
inline bool is_radian_between(const jereal &rad, const jereal &start_rad, const jereal &end_rad, const jereal &epsilon = Epsilon<jereal>()) noexcept
{
    auto sa = normalize_radian(start_rad);
    auto ea = normalize_radian(end_rad);
    auto r = normalize_radian(rad);
    if (fuzzy_equal(r, sa, epsilon) || fuzzy_equal(r, ea, epsilon)) {
        return true;
    }
    if (sa < ea) {
        return r > sa && r < ea;
    }
    else {
        return r > sa || r < ea;
    }
}

template<typename T>
inline bool is_point_equal(const T& x1, const T& y1, const T& x2, const T& y2, const T& epsilon = Epsilon<T>()) noexcept
{
    return fuzzy_equal(x1, x2, epsilon) && fuzzy_equal(y1, y2, epsilon);
}

/**
 * @brief 计算两个矢量的点积
 */
template<typename T>
inline jereal dot_product(const T &x1, const T &y1, const T &x2, const T &y2) noexcept
{
    return jemul(x1, x2) + jemul(y1, y2);
}

/**
 * @brief 计算两个矢量的叉积
 */
template<typename T>
inline jereal cross_product(const T &x1, const T &y1, const T &x2, const T &y2) noexcept
{
    return jemul(x1, y2) - jemul(y1, x2);
}

/**
 * @brief 矢量0->1, 0->2的叉积
 */
template<typename T>
inline jereal cross_product(const T &x0, const T &y0, const T &x1, const T &y1, const T &x2, const T &y2) noexcept
{
    return cross_product(x1 - x0, y1 - y0, x2 - x0, y2 - y0);
}

/**
 * @brief 两点平方距离
 * @param x1,y1: 点1坐标
 * @param x2,y2: 点2坐标
 * @return
 */
template<typename T>
inline jereal sqrdist(const T &x1, const T &y1, const T &x2, const T &y2) noexcept
{
    return jesqr(x2-x1) + jesqr(y2-y1);
}

/**
 * @brief 两点距离
 * @param x1,y1: 点1坐标
 * @param x2,y2: 点2坐标
 * @return
 */
template<typename T>
inline T distance(const T &x1, const T &y1, const T &x2, const T &y2) noexcept
{
    return jehypot<T>(x2 - x1, y2 - y1);
}

/**
 * @brief 从原点到点的弧度
 * @param x,y: 点坐标
 * @return
 */
template<typename T>
inline jereal point_radian(const T &x, const T &y) noexcept
{
    return std::atan2(y, x);
}

/**
 * @brief 从点1到点2构成的向量的弧度
 * @param x1,y1: 点1坐标
 * @param x2,y2: 点2坐标
 * @return
 */
template<typename T>
inline jereal vector_radian(const T &x1, const T &y1, const T &x2, const T &y2) noexcept
{
    return std::atan2(y2 - y1, x2 - x1);
}

/**
 * @brief 从起点位置延长多少倍后的点
 * @param x,y: 点坐标
 * @return
 */
template <typename T>
inline void project_point_ratio(const T&  xs, const T&  ys, const T& xe, const T& ye,
                                const jereal& t, T& nx, T& ny)
{
    nx = roundcast<T>(xs + jemul(t, (xe - xs)));
    ny = roundcast<T>(ys + jemul(t, (ye - ys)));
}

/**
 * @brief 从起点位置延长多少距离后的点
 * @param x,y: 点坐标
 * @return
 */
template <typename T>
inline void project_point_len(const T&  xs, const T&  ys, const T& xe, const T& ye,
                              const T& len, T& nx, T& ny)
{
    const T d = distance<T>(xs, ys, xe, ye);
    if (fuzzy_iszero(d)) {
        nx = xe;
        ny = ye;
        return;
    }
    jereal t = jediv(len, d);
    project_point_ratio<T>(xs, ys, xe, ye, t, nx, ny);
}

/**
 * @brief 判断三点是否共线
 * @param x1, y1 : 点1坐标
 * @param x2, y2 : 点2坐标
 * @param x3, y3 : 点3坐标
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool collinear(const T& x1, const T& y1, const T& x2, const T& y2, const T& x3, const T& y3, const T& epsilon = Epsilon<T>()) noexcept
{
    return fuzzy_equal<jereal>(jemul(x2 - x1, y3 - y1), jemul(x3 - x1, y2 - y1));
}


/**
 * @brief 判断两条线段的位置关系(v2位于v1的哪一侧或相交)
 * @param v1_xs, v1_ys, v1_xe, v1_ye : 线段1的起点和终点
 * @param v2_xs, v2_ys, v2_xe, v2_ye : 线段2的起点和终点
 * @param epsilon : 公差
 * @return 0:表示两线段相交; 1:表示v2在v1的左侧; -1:表示v2在v1的右侧
 */
template<typename T>
inline int orientation_vector2vector(
    const T &v1_xs, const T &v1_ys, const T &v1_xe, const T &v1_ye,
    const T &v2_xs, const T &v2_ys, const T &v2_xe, const T &v2_ye,
    const T &epsilon = Epsilon<T>()) noexcept
{
    // 判断两个向量的位置关系, 0表示两向量相交, 1表示v2在v1的左侧, -1表示v2在v1的右侧
    T a1 = v1_xe - v1_xs;
    T b1 = v1_ye - v1_ys;
    T a2 = v2_xe - v2_xs;
    T b2 = v2_ye - v2_ys;
    if (epsilon != 0.0) {
        if (fuzzy_iszero(a1, epsilon)) a1 = 0;
        if (fuzzy_iszero(b1, epsilon)) b1 = 0;
        if (fuzzy_iszero(a2, epsilon)) a2 = 0;
        if (fuzzy_iszero(b2, epsilon)) b2 = 0;
    }

    int sign1 = 0;
    int sign2 = 0;
    if (a1 && b2) {
        sign1 = ((a1 > 0) ^ (b2 > 0)) ? -1 : 1;
    }
    if (a2 && b1) {
        sign2 = ((a2 > 0) ^ (b1 > 0)) ? -1 : 1;
    }

    if (sign1 != sign2) {
        return (sign1 < sign2) ? -1 : 1;
    }

    jereal a3 = jemul(a1 < 0 ? -a1 : a1, b2 < 0 ? -b2 : b2);
    jereal b3 = jemul(b1 < 0 ? -b1 : b1, a2 < 0 ? -a2 : a2);
    if (fuzzy_equal<jereal>(a3, b3, epsilon)) {
        return 0;
    }
    return ((a3 < b3) ^ (sign1 == 1)) ? 1 : -1;
}

/**
 * @brief 判断点与直线的位置关系(点位于v1的哪一侧,或点在直线上)
 * @param vx1, vy1, vx2, vy2 : 线段的两个端点
 * @param px, py : 点的坐标
 * @param epsilon : 公差
 * @return 0:表示点线段及其两端延长线上; 1:表示点在线段左侧; -1:表示点在线段右侧
 */
template<typename T>
inline int orientation_vector2point(const T& vx1, const T& vy1, const T& vx2, const T& vy2,
                                    const T& px, const T& py,
                                    const T& epsilon = Epsilon<T>()) noexcept
{
    return orientation_vector2vector<T>(vx1, vy1, vx2, vy2, vx2, vy2, px, py, epsilon);
}

/**
 * @brief 点与圆的位置关系
 * @param cx, cy, cr : 圆心和半径
 * @param px, py : 点的坐标
 * @param epsilon : 公差
 * @return 0:表示点在圆上; 1:表示点在圆内; -1:表示点在圆外
*/
template<typename T>
inline int orientation_circle2point(const T& cx, const T& cy, const T&cr,
                                    const T& px, const T& py,
                                    const T& epsilon = Epsilon<T>()) noexcept
{

    T d = distance<T>(cx, cy, px, py);
    if (fuzzy_equal<T>(d, cr, epsilon)) {
        return 0;
    }
    return (d < cr) ? 1 : -1;
}


/**
 * @brief 判断点与矩形的位置关系
 * @param rect_x1, rect_y1, rect_x2, rect_y2 : 矩形的对角线两个点
 * @param px, py : 点的坐标
 * @param epsilon : 公差
 * @return 0:表示点在矩形边界上; 1:表示点在矩形内; -1:表示点在矩形外
 */
template<typename T>
inline int orientation_rect2point(const T& rect_x1, const T& rect_y1, const T& rect_x2, const T& rect_y2,
                                  const T& px, const T& py,
                                  const T& epsilon = Epsilon<T>()) noexcept
{
    T x1 = jemin(rect_x1, rect_x2);
    T x2 = jemax(rect_x1, rect_x2);
    T y1 = jemin(rect_y1, rect_y2);
    T y2 = jemax(rect_y1, rect_y2);
    //判断是否在矩形边界上
    if (((fuzzy_lessequal(x1, px, epsilon) && fuzzy_lessequal(px, x2, epsilon)) && ((fuzzy_equal(py, y1, epsilon) || fuzzy_equal(py, y2, epsilon))))
        || ((fuzzy_lessequal(y1, py, epsilon) && fuzzy_lessequal(py, y2, epsilon)) && ((fuzzy_equal(px, x1, epsilon) || fuzzy_equal(px, x2, epsilon)))))
    {
        return 0;
    }
    return (x1 < px && px < x2 && y1 < py && py < y2) ? 1 : -1;
}


/**
 * @brief 通过凸度计算圆弧的圆心
 * @param circle_x, circle_y, circle_r : 圆心和半径
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段起点和终点
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T>
inline void arc_center_from_bulge(const T &xs, const T &ys, const T &xe, const T &ye, const jereal &bulge,
                                  T &center_x, T &center_y ) noexcept
{
    jereal b = 0.5 * (1.0 / bulge - bulge);
    center_x = roundcast<T, jereal>( 0.5 * ((xs + xe) - b * (ye - ys)));
    center_y = roundcast<T, jereal>( 0.5 * ((ys + ye) + b * (xe - xs)));
}

/**
 * @brief 通过圆弧角度计算圆弧的圆心
 * @param xs, ys, xe, ye : 线段起点和终点
 * @param angle : 圆弧角度(正为逆时针,负为顺时针)
 * @param center_x, center_y : 返回的圆心坐标
 */
template<typename T>
inline void arc_center_from_angle(const T &xs, const T &ys, const T &xe, const T &ye, const jereal &angle,
                                  T &center_x, T &center_y) noexcept
{
    jereal bulge = deg2bulge(angle);
    arc_center_from_bulge<T>(xs, ys, xe, ye, bulge, center_x, center_y);
}

/**
 * @brief 计算圆弧的弧度
 * @param xs, ys, xe, ye ,xc, yc, cw : 圆弧的起点,终点,圆心,是否顺时针
 * @return 圆弧的夹角弧度
 */
template<typename T>
inline jereal arc_radian_from_secw(const T& xs, const T& ys, const T& xe, const T& ye, const T& xc, const T& yc, bool cw) noexcept
{
    jereal sa = std::atan2(yc - ys, xs - xc);
    jereal ea = std::atan2(yc - ye, xe - xc);
    if (cw) {
        if (ea <= sa) ea += PI_MUL_2;
    }
    else {
        if (sa <= ea) sa += PI_MUL_2;
    }
    return sa - ea;
}

/**
 * @brief 计算圆弧的角度
 * @param xs, ys, xe, ye ,xc, yc, cw : 圆弧的起点,终点,圆心,是否顺时针
 * @return 圆弧的夹角角度
 */
template<typename T>
inline jereal arc_degree_from_secw(const T& xs, const T& ys, const T& xe, const T& ye, const T& xc, const T& yc, bool cw) noexcept
{
    return rad2deg(arc_radian_from_secw(xs, ys, xe, ye, xc, yc, cw));
}

/**
 * @brief 判断点是否在直线上
 * @param px, py : 点坐标
 * @param line_x1, line_y1, line_x2, line_y2: 直线的两点
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_line(const T& px, const T& py,
                             const T& line_x1, const T& line_y1, const T& line_x2, const T& line_y2,
                             const T& epsilon = Epsilon<T>()) noexcept
{
    return collinear<T>(line_x1, line_y1, line_x2, line_y2, px, py, epsilon);
}

/**
 * @brief 点到直线的最近点
 * @param px, py : 点坐标
 * @param line_x1, line_y1, line_x2, line_y2 : 直线起点和终点
 * @param closest_x, closest_y : 返回的最近点坐标
 * @return 如果点在直线上及其两端的延长线上,则没有最近点, 返回0, 否则返回1
 */
template<typename T>
inline void closest_point2line(const T &px, const T &py,
                               const T &line_x1, const T &line_y1, const T &line_x2, const T &line_y2,
                               T &closest_x, T &closest_y) noexcept
{
    const T vx = line_x2 - line_x1;
    const T vy = line_y2 - line_y1;
    const T wx = px - line_x1;
    const T wy = py - line_y1;
    const jereal c1 = jemul(vx, wx) + jemul(vy, wy);
    const jereal c2 = jemul(vx, vx) + jemul(vy, vy);
    const jereal t = c1 / c2;
    closest_x = roundcast<T, jereal>(line_x1 + t * vx);
    closest_y = roundcast<T, jereal>(line_y1 + t * vy);
}


/**
 * @brief 判断点是否在线段上
 * @param px, py : 点坐标
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段起点和终点
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_segment(const T& px, const T& py,
                                const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                const T& epsilon = Epsilon<T>()) noexcept
{
    if (((fuzzy_lessequal<T>(seg_xs, px, epsilon) && fuzzy_lessequal<T>(px, seg_xe, epsilon))
         || (fuzzy_lessequal<T>(seg_xe, px, epsilon) && fuzzy_lessequal<T>(px, seg_xs, epsilon)))
        && ((fuzzy_lessequal<T>(seg_ys, py, epsilon) && fuzzy_lessequal<T>(py, seg_ye, epsilon))
            || (fuzzy_lessequal<T>(seg_ye, py, epsilon) && fuzzy_lessequal<T>(py, seg_ys, epsilon))))
    {
        return collinear<T>(seg_xs, seg_ys, seg_xe, seg_ye, px, py, epsilon);
    }
    return false;
}

/**
 * @brief 判断点是否在射线上
 * @param px, py : 点坐标
 * @param ray_xo, ray_yo : 射线起点
 * @param ray_xd, ray_yd : 射线方向
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_ray(const T& px, const T& py,
                            const T& ray_xo, const T& ray_yo, const T& ray_xd, const T& ray_yd,
                            const T& epsilon = Epsilon<T>()) noexcept
{
    const jereal t = dot_product<jereal>(ray_xd, ray_yd, px - ray_xo, py - ray_yo);
    if (fuzzy_greaterequal<jereal>(t, 0, epsilon)) {
        T nx = roundcast<T, jereal>(ray_xo + t * ray_xd);
        T ny = roundcast<T, jereal>(ray_yo + t * ray_yd);
        return fuzzy_equal<T>(px, nx, epsilon) && fuzzy_equal<T>(py, ny, epsilon);
    }
    return false;
}


/**
 * @brief 判断点是否在矩形上
 * @param px, py : 点坐标
 * @param rect_x1, rect_y1, rect_x2, rect_y2 : 矩形的对角线两个点
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_rect(const T& px, const T& py,
                             const T& rect_x1, const T& rect_y1, const T& rect_x2, const T& rect_y2,
                             const T& epsilon = Epsilon<T>()) noexcept
{
    T x1 = jemin(rect_x1, rect_x2);
    T x2 = jemax(rect_x1, rect_x2);
    T y1 = jemin(rect_y1, rect_y2);
    T y2 = jemax(rect_y1, rect_y2);
    return ((fuzzy_lessequal(x1, px, epsilon) && fuzzy_lessequal(px, x2, epsilon)) && (fuzzy_equal(py, y1, epsilon) || fuzzy_equal(py, y2, epsilon)))
           || ((fuzzy_lessequal(y1, py, epsilon) && fuzzy_lessequal(py, y2, epsilon)) && (fuzzy_equal(px, x1, epsilon) || fuzzy_equal(px, x2, epsilon)));
}

/**
 * @brief 判断点是否在圆上
 * @param px, py : 点坐标
 * @param cx, cy : 圆心坐标
 * @param r : 圆半径
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_circle(const T& px, const T& py,
                               const T& cx, const T& cy, const T& r,
                               const T& epsilon = Epsilon<T>()) noexcept
{
    return fuzzy_equal<jereal>(sqrdist<jereal>(px, py, cx, cy), jesqr(r), epsilon);
}


/**
 * @brief 判断点是否在圆弧上(圆弧定义为圆心(c),半径(r),起始弧度,结束弧度(a))
 * @param px, py : 点坐标
 * @param arc_cx, arc_cy, arc_r : 圆弧的圆心和半径
 * @param arc_start_rad, arc_end_rad : 圆弧的起始和结束弧度
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool is_point_on_arc_cra(const T& px, const T& py,
                                const T &arc_cx, const T &arc_cy, const T &arc_r, const T &arc_start_rad, const T &arc_end_rad,
                                const T& epsilon = Epsilon<T>()) noexcept
{
    if (!is_point_on_circle<T>(px, py, arc_cx, arc_cy, arc_r, epsilon)) {
        return false;
    }
    jereal point_rad = atan2(py - arc_cy, px - arc_cx);
    return is_radian_between(point_rad, arc_start_rad, arc_end_rad);
}

/**
 * @brief 判断点是否在圆弧上(圆弧定义为起点(s),终点(e),圆心(c),是否顺时针(w))
 * @param px, py : 点坐标
 * @param arc_xs, arc_ys, arc_xe, arc_ye : 圆弧的起点和终点
 * @param arc_xc, arc_yc : 圆弧的圆心
 * @param arc_cw : 圆弧是否顺时针
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_arc_secw(const T &px, const T &py,
                                 const T &arc_xs, const T &arc_ys, const T &arc_xe, const T &arc_ye,
                                 const T &arc_xc, const T &arc_yc, bool arc_cw,
                                 const T& epsilon = Epsilon<T>())
{
    T arc_r = jehypot(arc_xc - arc_xs, arc_yc - arc_ys);
    jereal sa = vector_radian<T>(arc_xc, arc_yc, arc_xs, arc_ys);
    jereal ea = vector_radian<T>(arc_xc, arc_yc, arc_xe, arc_ye);
    if (arc_cw) {
        std::swap(sa, ea);
    }
    return is_point_on_arc_cra<T>(px, py, arc_xc, arc_yc, arc_r, sa, ea, epsilon);
}


/**
 * @brief 判断点是否在圆弧上(圆弧定义为起点(s),终点(e),角度(a))
 * @param px, py : 点坐标
 * @param arc_xs, arc_ys, arc_xe, arc_ye : 圆弧的起点和终点
 * @param arc_angle : 圆弧的角度
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_point_on_arc_sea(const T &px, const T &py,
                                const T &arc_xs, const T &arc_ys, const T &arc_xe, const T &arc_ye, const jereal &arc_angle
                                )
{
    T arc_xc, arc_yc;
    arc_center_from_angle<T>(arc_xs, arc_ys, arc_xe, arc_ye, arc_angle, arc_xc, arc_yc);
    bool arc_cw = arc_angle < 0;
    T arc_r = jehypot(arc_xc - arc_xs, arc_yc - arc_ys);
    jereal sa = vector_radian<T>(arc_xc, arc_yc, arc_xs, arc_ys);
    jereal ea = vector_radian<T>(arc_xc, arc_yc, arc_xe, arc_ye);
    if (arc_cw) {
        std::swap(sa, ea);
    }
    return is_point_on_arc_cra<T>(px, py, arc_xc, arc_yc, arc_r, sa, ea);
}


/**
 * @brief 判断点是否在圆内
 * @param px, py : 点坐标
 * @param cx, cy, cr : 圆心和半径
 * @param consider_touch : 是否考虑只要碰到就算在圆内
 * @param epsilon : 公差
 */
template<typename T>
inline bool contains_point_in_circle(const T& px, const T& py,
                                   const T& cx, const T& cy, const T& cr,
                                   const bool &consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    int ori = orientation_circle2point<T>(cx, cy, cr, px, py, epsilon);
    if (consider_touch) {
        return ori >= 0;
    }
    else {
        return ori > 0;
    }
}

/**
 * @brief 判断点是否在矩形内
 * @param px, py : 点坐标
 * @param rect_x1, rect_y1, rect_x2, rect_y2 : 矩形的对角线两个点
 * @param consider_touch : 是否考虑只要碰到就算在矩形内
 * @param epsilon : 公差
 */
template<typename T>
inline bool contains_point_in_rect(const T& px, const T& py,
                                 const T& rect_x1, const T& rect_y1, const T& rect_x2, const T& rect_y2,
                                 const bool &consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    int ori = orientation_rect2point<T>(rect_x1, rect_y1, rect_x2, rect_y2, px, py, epsilon);
    if (consider_touch) {
        return ori >= 0;
    }
    else {
        return ori > 0;
    }
}

template<typename T>
inline bool contains_point_in_monrect(const T& px, const T& py,
                                        const T& rect_xmin, const T& rect_ymin, const T& rect_xmax, const T& rect_ymax,
                                        const bool &consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (consider_touch) {
        return (fuzzy_lessequal(rect_xmin, px, epsilon) && fuzzy_lessequal(px, rect_xmax, epsilon)
                && fuzzy_lessequal(rect_ymin, py, epsilon) && fuzzy_lessequal(py, rect_ymax, epsilon));
    }
    return (rect_xmin < px && px < rect_xmax && rect_ymin < py && py < rect_ymax);
}

template<typename T>
inline bool contains_monrect_in_monrect(const T& recta_xmin, const T& recta_ymin, const T& recta_xmax, const T& recta_ymax,
                                        const T& rectb_xmin, const T& rectb_ymin, const T& rectb_xmax, const T& rectb_ymax,
                                        const bool &consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    // 判断矩形A是否完全包含在矩形B内
    if (consider_touch) {
        return (fuzzy_lessequal(rectb_xmin, recta_xmin, epsilon) && fuzzy_lessequal(recta_xmax, rectb_xmax, epsilon)
                && fuzzy_lessequal(rectb_ymin, recta_ymin, epsilon) && fuzzy_lessequal(recta_ymax, rectb_ymax, epsilon));
    }
    return (rectb_xmin < recta_xmin && recta_xmax < rectb_xmax && rectb_ymin < recta_ymin && recta_ymax < rectb_ymax);
}

/**
 * @brief 判断 a 矩形是否被 b 矩形包含
 */
template<typename T>
inline bool contains_rect_in_rect(const T& recta_x1, const T& recta_y1, const T& recta_x2, const T& recta_y2,
                               const T& rectb_x1, const T& rectb_y1, const T& rectb_x2, const T& rectb_y2,
                               bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    T aminx = jemin(recta_x1, recta_x2);
    T amaxx = jemax(recta_x1, recta_x2);
    T aminy = jemin(recta_y1, recta_y2);
    T amaxy = jemax(recta_y1, recta_y2);
    T bminx = jemin(rectb_x1, rectb_x2);
    T bmaxx = jemax(rectb_x1, rectb_x2);
    T bminy = jemin(rectb_y1, rectb_y2);
    T bmaxy = jemax(rectb_y1, rectb_y2);

    if (consider_touch && fuzzy_lessequal(bminx, aminx, epsilon) && fuzzy_greaterequal(bmaxx, amaxx, epsilon)
            && fuzzy_lessequal(bminy, aminy, epsilon) && fuzzy_greaterequal(bmaxy, amaxy, epsilon)) {
        return true;
    } else if (!consider_touch && bminx < aminx && bmaxx > amaxx && bminy < aminy && bmaxy > amaxy) {
        return true;
    }
    return false;
}


/**
 * @brief 判断值是否在区间内
 * @param v : 待判断的值
 * @param minv : 区间最小值
 * @param maxv : 区间最大值
 * @param consider_touch : 是否考虑只要碰到就算在区间内
 * @param epsilon : 公差
 */
template<typename T>
inline bool is_value_between(const T& v, const T& minv, const T& maxv, bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (consider_touch) {
        return fuzzy_greaterequal<T>(v, minv, epsilon) && fuzzy_lessequal<T>(v, maxv, epsilon);
    }
    return (v > minv) && (v < maxv);
}

template<typename T>
inline bool is_interval_cross(const T& amin, const T& amax, const T& bmin, const T& bmax, bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (consider_touch) {
        return fuzzy_lessequal<T>(amin, bmax, epsilon) && fuzzy_greaterequal<T>(amax, bmin, epsilon);
    }
    return (amin < bmax) && (amax > bmin);
}


/**
 * @brief 判断两个区间是否相交
 * @param a1, a2 : 区间1的两个值
 * @param b1, b2 : 区间2的两个值
 * @param consider_touch : 是否考虑只要碰到就算相交
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool intersects_interval2interval(const T &a1, const T& a2, const T& b1, const T& b2,
                                         bool consider_touch = true, const T &epsilon = Epsilon<T>()) noexcept
{
    return is_interval_cross<T>(jemin(a1, a2), jemax(a1, a2), jemin(b1, b2), jemax(b1, b2), consider_touch, epsilon);
}

/**
 * @brief 判断两个矩形是否相交
 * @param recta_x1, recta_y1, recta_x2, recta_y2 : 矩形1的两个对角线坐标
 * @param rectb_x1, rectb_y1, rectb_x2, rectb_y2 : 矩形2的两个对角线坐标
 * @param consider_touch : 是否考虑只要碰到就算相交
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool intersects_rect2rect(const T& recta_x1, const T& recta_y1, const T& recta_x2, const T& recta_y2,
                                 const T& rectb_x1, const T& rectb_y1, const T& rectb_x2, const T& rectb_y2,
                                 bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    return intersects_interval2interval<T>(recta_x1, recta_x2, rectb_x1, rectb_x2, consider_touch, epsilon) &&
           intersects_interval2interval<T>(recta_y1, recta_y2, rectb_y1, rectb_y2, consider_touch, epsilon);
}


/**
 * @brief 判断两个直线是否相交
 * @param linea_x1, linea_y1, linea_x2, linea_y2 : 直线1的两个点
 * @param lineb_x1, lineb_y1, lineb_x2, lineb_y2 : 直线2的两个点
 * @param consider_touch : 是否考虑只要碰到就算相交
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool intersects_line2line(const T& linea_x1, const T& linea_y1, const T& linea_x2, const T& linea_y2,
                                 const T& lineb_x1, const T& lineb_y1, const T& lineb_x2, const T& lineb_y2,
                                 bool consider_touch = false,
                                 const T& epsilon = Epsilon<T>()) noexcept
{
    if (fuzzy_notequal(jemul(linea_x1 - linea_x2, lineb_y1 - lineb_y2), jemul(linea_y1 - linea_y2, lineb_x1 - lineb_x2), epsilon)) {
        return true;
    }
    else if (consider_touch && collinear<T>(linea_x1, linea_y1, linea_x2, linea_y2, lineb_x1, lineb_y1, epsilon)) {
        return true;
    }
    return false;
}



/**
 * @brief 判断直线与线段是否相交
 * @param line_x1, line_y1, line_x2, line_y2 : 直线的两个点
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段的起点和终点
 * @param consider_touch : 是否考虑只要碰到就算相交
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool intersects_line2segment(const T& line_x1, const T& line_y1, const T& line_x2, const T& line_y2,
                                    const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                    bool consider_touch = false,
                                    const T& epsilon = Epsilon<T>()) noexcept
{
    int ora_1 = orientation_vector2point<T>(line_x1, line_y1, line_x2, line_y2, seg_xs, seg_ys, epsilon);
    int ora_2 = orientation_vector2point<T>(line_x1, line_y1, line_x2, line_y2, seg_xe, seg_ye, epsilon);
    if (ora_1 * ora_2 > 0) {
        return false;
    }

    if (consider_touch || (ora_1 && ora_2)) {
        return true;
    }
    return false;
}

/**
 * @brief 判断两个线段是否相交
 * @param sega_xs, sega_ys, sega_xe, sega_ye : 线段1的起点和终点
 * @param segb_xs, segb_ys, segb_xe, segb_ye : 线段2的起点和终点
 * @param consider_touch : 是否考虑只要碰到就算相交(注意:这里的touch为false时,只会排除线段首尾想接的情况,
 *                         如果是在一个线段的端点位于另一个线段上(但不是端点位置),不会被排除)
 * @param epsilon : 公差
 * @return
 */
template<typename T>
inline bool intersects_segment2segment(const T& sega_xs, const T& sega_ys, const T& sega_xe, const T& sega_ye,
                                       const T& segb_xs, const T& segb_ys, const T& segb_xe, const T& segb_ye,
                                       bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (!intersects_rect2rect<T>(sega_xs, sega_ys, sega_xe, sega_ye, segb_xs, segb_ys, segb_xe, segb_ye, true, epsilon)) {
        return false;
    }
    int ora_1 = orientation_vector2point<T>(sega_xs, sega_ys, sega_xe, sega_ye, segb_xs, segb_ys, epsilon);
    int ora_2 = orientation_vector2point<T>(sega_xs, sega_ys, sega_xe, sega_ye, segb_xe, segb_ye, epsilon);
    if (ora_1 * ora_2 > 0) {
        return false;
    }

    int orb_1 = orientation_vector2point<T>(segb_xs, segb_ys, segb_xe, segb_ye, sega_xs, sega_ys, epsilon);
    int orb_2 = orientation_vector2point<T>(segb_xs, segb_ys, segb_xe, segb_ye, sega_xe, sega_ye, epsilon);
    if (orb_1 * orb_2 > 0) {
        return false;
    }

    if (consider_touch || (ora_1 && ora_2) || (orb_1 && orb_2)) {
        return true;
    }

    if (ora_1 || ora_2) {
        return false;
    }

    return intersects_interval2interval<T>(sega_xs, sega_xe, segb_xs, segb_xe, consider_touch, epsilon) &&
           intersects_interval2interval<T>(sega_ys, sega_ye, segb_ys, segb_ye, consider_touch, epsilon);

}

template<typename T>
inline bool intersects_segment2rect(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                    const T& rect_x1, const T& rect_y1, const T& rect_x2, const T& rect_y2,
                                    bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (!intersects_rect2rect<T>(seg_xs, seg_ys, seg_xe, seg_ye, rect_x1, rect_y1, rect_x2, rect_y2, true, epsilon)) {
        return false;
    }

    if (intersects_segment2segment<T>(seg_xs, seg_ys, seg_xe, seg_ye, rect_x1, rect_y1, rect_x1, rect_y2, consider_touch, epsilon)) {
        return true;
    }
    if (intersects_segment2segment<T>(seg_xs, seg_ys, seg_xe, seg_ye, rect_x1, rect_y1, rect_x2, rect_y1, consider_touch, epsilon)) {
        return true;
    }
    if (intersects_segment2segment<T>(seg_xs, seg_ys, seg_xe, seg_ye, rect_x2, rect_y2, rect_x2, rect_y1, consider_touch, epsilon)) {
        return true;
    }
    if (intersects_segment2segment<T>(seg_xs, seg_ys, seg_xe, seg_ye, rect_x2, rect_y2, rect_x1, rect_y2, consider_touch, epsilon)) {
        return true;
    }
    return false;
}

template<typename T>
inline bool intersects_circle2circle(const T& circle1_x, const T& circle1_y, const T& circle1_r,
                                     const T& circle2_x, const T& circle2_y, const T& circle2_r,
                                     bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (consider_touch) {
        return fuzzy_lessequal<T>(sqrdist<T>(circle1_x, circle1_y, circle2_x, circle2_y), jesqr(circle1_r + circle2_r), epsilon);
    }
    return sqrdist<T>(circle1_x, circle1_y, circle2_x, circle2_y) < jesqr(circle1_r + circle2_r);
}



template<typename T>
inline bool intersects_line2circle(const T& line_x1, const T& line_y1, const T& line_x2, const T& line_y2,
                                   const T& circle_x, const T& circle_y, const T& circle_r,
                                   bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{

    T closest_x, closest_y;
    closest_point2line<T>(circle_x, circle_y, line_x1, line_y1, line_x2, line_y2, closest_x, closest_y);
    return contains_point_in_circle<T>(closest_x, closest_y, circle_x, circle_y, circle_r, consider_touch, epsilon);
}

template<typename T>
inline bool intersects_segment2circle(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                      const T& circle_x, const T& circle_y, const T& circle_r,
                                      bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    int ori1 = orientation_circle2point<T>(circle_x, circle_y, circle_r, seg_xs, seg_ys, epsilon);
    int ori2 = orientation_circle2point<T>(circle_x, circle_y, circle_r, seg_xe, seg_ye, epsilon);

    if (consider_touch && (ori1 == 0 || ori2 == 0)) {
        return true;
    }
    // 一个点在圆内,一个点在圆外, 一定相交
    if (ori1 * ori2 < 0) {
        return true;
    }
    // 两个点都在圆内, 一定不相交
    if (ori1 > 0 && ori2 > 0) {
        return false;
    }
    // 两个点都在圆外, 则判断圆心到直线的距离是否小于半径
    return intersects_line2circle<T>(seg_xs, seg_ys, seg_xe, seg_ye, circle_x, circle_y, circle_r, consider_touch, epsilon);
}

template<typename T>
inline int intersection_line2line(const T& line1_x1, const T& line1_y1, const T& line1_x2, const T& line1_y2,
                                  const T& line2_x1, const T& line2_y1, const T& line2_x2, const T& line2_y2,
                                  T& ix, T& iy)
{
    T dx1 = line1_x2 - line1_x1;
    T dy1 = line1_y2 - line1_y1;
    T dx2 = line2_x2 - line2_x1;
    T dy2 = line2_y2 - line2_y1;
    jereal d = cross_product(dx1, dy1, dx2, dy2);
    if (fuzzy_equal<jereal>(d, 0)) {
        // 平行或共线
        return 0;
    }

    T dx3 = line2_x1 - line1_x1;
    T dy3 = line2_y1 - line1_y1;
    jereal t = cross_product(dx3, dy3, dx2, dy2) / d;
    ix = roundcast<T, jereal>(line1_x1 + t * dx1);
    iy = roundcast<T, jereal>(line1_y1 + t * dy1);
    return 1;
}

template<typename T>
inline int intersection_segment2segment(const T& seg1_xs, const T& seg1_ys, const T& seg1_xe, const T& seg1_ye,
                                         const T& seg2_xs, const T& seg2_ys, const T& seg2_xe, const T& seg2_ye,
                                         T& ix, T& iy)
{
    T dx1 = seg1_xe - seg1_xs;
    T dy1 = seg1_ye - seg1_ys;
    T dx2 = seg2_xe - seg2_xs;
    T dy2 = seg2_ye - seg2_ys;
    jereal d = cross_product(dx1, dy1, dx2, dy2);
    if (fuzzy_equal<jereal>(d, 0)) {
        // 平行或共线
        return -1;
    }

    T dx3 = seg2_xs - seg1_xs;
    T dy3 = seg2_ys - seg1_ys;
    jereal t1 = cross_product(dx3, dy3, dx2, dy2) / d;
    jereal t2 = cross_product(dx3, dy3, dx1, dy1) / d;
    if (t1 >= 0 && t1 <= 1 && t2 >= 0 && t2 <= 1) {
        ix = roundcast<T, jereal>(seg1_xs + t1 * dx1);
        iy = roundcast<T, jereal>(seg1_ys + t1 * dy1);
        return 1;
    }
    // 不相交
    return 0;
}



template<typename T>
inline int intersection_line2circle(const T& line_x1, const T& line_y1, const T& line_x2, const T& line_y2,
                                    const T& circle_x, const T& circle_y, const T& circle_r,
                                    T& inter_x1, T& inter_y1, T& inter_x2, T& inter_y2,
                                    const T& epsilon = Epsilon<T>()) noexcept
{
    const jereal a = jesqr(line_x2 - line_x1) + jesqr(line_y2 - line_y1);
    const jereal b = 2.0 * (jemul(line_x2 - line_x1, line_x1 - circle_x) + jemul(line_y2 - line_y1, line_y1 - circle_y));
    const jereal c = jesqr(circle_x) + jesqr(circle_y) + jesqr(line_x1) + jesqr(line_y1) - 2.0 * (jemul(circle_x, line_x1) + jemul(circle_y, line_y1)) - jesqr(circle_r);
    const jereal delta = jesqr(b) - 4.0 * a * c;
    if (delta < 0) {
        return 0;
    }
    else if (fuzzy_iszero<jereal>(delta, epsilon)) {
        jereal t = -0.5 * b / a;
        inter_x1 = roundcast<T, jereal>(line_x1 + t * (line_x2 - line_x1));
        inter_y1 = roundcast<T, jereal>(line_y1 + t * (line_y2 - line_y1));
        return 1;
    }
    else {
        const jereal sqrt_delta = std::sqrt(delta);
        jereal t1 = (-b + sqrt_delta) / (2.0 * a);
        jereal t2 = (-b - sqrt_delta) / (2.0 * a);
        inter_x1 = roundcast<T, jereal>(line_x1 + t1 * (line_x2 - line_x1));
        inter_y1 = roundcast<T, jereal>(line_y1 + t1 * (line_y2 - line_y1));
        inter_x2 = roundcast<T, jereal>(line_x1 + t2 * (line_x2 - line_x1));
        inter_y2 = roundcast<T, jereal>(line_y1 + t2 * (line_y2 - line_y1));
        return 2;
    }
}


template<typename T>
inline int intersection_segment2circle(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                       const T& cx, const T& cy, const T& cr,
                                       T& inter_x1, T& inter_y1, T& inter_x2, T& inter_y2,
                                       const T& epsilon = Epsilon<T>()) noexcept
{
    T ix1, iy1, ix2, iy2;
    int inter_count = intersection_line2circle(seg_xs, seg_ys, seg_xe, seg_ye, cx, cy, cr, ix1, iy1, ix2, iy2, epsilon);
    if (inter_count <= 0) return 0;

    int res_count = 0;
    if (inter_count >= 1) {
        if (is_point_on_segment(ix1, iy1, seg_xs, seg_ys, seg_xe, seg_ye, epsilon)) {
            inter_x1 = ix1;
            inter_y1 = iy1;
            res_count++;
        }
        if(inter_count == res_count)
            return -1;
    }
    if (inter_count >= 2) {
        if (is_point_on_segment(ix2, iy2, seg_xs, seg_ys, seg_xe, seg_ye, epsilon)) {
            if (res_count > 0) {
                inter_x2 = ix2;
                inter_y2 = iy2;
            }
            else {
                inter_x1 = ix2;
                inter_y1 = iy2;
            }
            res_count++;
        }
    }
    return res_count;
}


template<typename T>
inline int intersection_segment2arc(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                    const T& arc_xc, const T& arc_yc,
                                    const T& arc_radius, const T& arc_start_rad, const T& arc_end_rad,
                                    T& inter_x1, T& inter_y1, T& inter_x2, T& inter_y2,
                                    const T& epsilon = Epsilon<T>()) noexcept
{
    T ci_x1, ci_y1, ci_x2, ci_y2;
    int ci_inter_num = intersection_segment2circle<T>(seg_xs, seg_ys, seg_xe, seg_ye, arc_xc, arc_yc, arc_radius, ci_x1, ci_y1, ci_x2, ci_y2, epsilon);
    if (ci_inter_num == 0) return 0;
    int arc_inter_count = 0;
    if (ci_inter_num >= 1) {
        if (is_radian_between(vector_radian<T>(arc_xc, arc_yc, ci_x1, ci_y1), arc_start_rad, arc_end_rad, epsilon)) {
            inter_x1 = ci_x1;
            inter_y1 = ci_y1;
            arc_inter_count++;
        }
        if (ci_inter_num == 1 && arc_inter_count == 1) {
            return -1;
        }
    }
    if (ci_inter_num >= 2) {
        if (is_radian_between(vector_radian<T>(arc_xc, arc_yc, ci_x2, ci_y2), arc_start_rad, arc_end_rad, epsilon)) {
            if (arc_inter_count > 0) {
                inter_x2 = ci_x2;
                inter_y2 = ci_y2;
            }
            else {
                inter_x1 = ci_x2;
                inter_y1 = ci_y2;
            }
            arc_inter_count++;
        }
    }
    return arc_inter_count;
}


template<typename T>
inline int intersection_segment2monarc(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                       const T& arc_xs, const T& arc_ys, const T& arc_xe, const T& arc_ye,
                                       const T& arc_xc, const T& arc_yc, const T& arc_radius,
                                       T& inter_x1, T& inter_y1, T& inter_x2, T& inter_y2,
                                       const T& epsilon = Epsilon<T>()) noexcept
{

    T ci_x1, ci_y1, ci_x2, ci_y2;
    int ci_inter_num = intersection_segment2circle<T>(seg_xs, seg_ys, seg_xe, seg_ye, arc_xc, arc_yc, arc_radius, ci_x1, ci_y1, ci_x2, ci_y2, epsilon);
    if (ci_inter_num == 0) return 0;
    int arc_inter_count = 0;
    T aminx = jemin(arc_xs, arc_xe);
    T amaxx = jemax(arc_xs, arc_xe);
    T aminy = jemin(arc_ys, arc_ye);
    T amaxy = jemax(arc_ys, arc_ye);
    if (jeabs(ci_inter_num) >= 1) {
        if (is_value_between(ci_x1, aminx, amaxx, true, epsilon) && is_value_between(ci_y1, aminy, amaxy, true, epsilon)) {
            inter_x1 = ci_x1;
            inter_y1 = ci_y1;
            arc_inter_count++;
        }
        if (ci_inter_num < 2 && arc_inter_count == 1) {
            return ci_inter_num;
        }
    }
    if (ci_inter_num >= 2) {
        if (is_value_between(ci_x2, aminx, amaxx, true, epsilon) && is_value_between(ci_y2, aminy, amaxy, true, epsilon)) {
            if (arc_inter_count > 0) {
                inter_x2 = ci_x2;
                inter_y2 = ci_y2;
            }
            else {
                inter_x1 = ci_x2;
                inter_y1 = ci_y2;
            }
            arc_inter_count++;
        }
    }
    return arc_inter_count;
}

/**
 * @brief 获取水平射线和单调弧的交点
 */
template<typename T>
inline int intersection_rray2monarc(const T &ray_x, const T &ray_y,
                                    const T &arc_sx, const T &arc_sy, const T &arc_ex, const T &arc_ey,
                                    const T &arc_cx, const T &arc_cy, const T &arc_r,
                                    T& inter_x1, T& inter_y1,
                                    const T& epsilon = Epsilon<T>()) noexcept
{
    int arc_inter_count = 0;
    T aminx = jemin(arc_sx, arc_ex);
    T amaxx = jemax(arc_sx, arc_ex);
    T aminy = jemin(arc_sy, arc_ey);
    T amaxy = jemax(arc_sy, arc_ey);


    if (ray_x > amaxx || ray_y < aminy || ray_y > amaxy) {
        return arc_inter_count;
    }

    T inter_minx = jemax(aminx, ray_x);
    T dx = jesqrt<T>(jesqr(arc_r) - jesqr(ray_y - arc_cy));
    T arc_x1 = arc_cx + dx;
    T arc_x2 = arc_cx - dx;
    if (is_value_between(arc_x1, inter_minx, amaxx, epsilon)) {
        inter_x1 = arc_x1;
        inter_y1 = ray_y;
        arc_inter_count++;
    } else if (is_value_between(arc_x2, inter_minx, amaxx, epsilon)) {
        inter_x1 = arc_x2;
        inter_y1 = ray_y;
        arc_inter_count++;
    }

    return arc_inter_count;
}


template<typename T>
inline int intersection_circle2circle(const T& circle1_x, const T& circle1_y, const T& circle1_r,
                                      const T& circle2_x, const T& circle2_y, const T& circle2_r,
                                      T& inter_x1, T& inter_y1, T& inter_x2, T& inter_y2,
                                      const T& epsilon = Epsilon<T>()) noexcept
{
    T dx = circle2_x - circle1_x;
    T dy = circle2_y - circle1_y;
    jereal d = jehypot<jereal>(dx, dy);
    if (d > circle1_r + circle2_r || d < jeabs(circle1_r - circle2_r)) {
        // 圆不相交
        return 0;
    }
    else if (fuzzy_iszero(roundcast<T>(d), epsilon) && fuzzy_equal(circle1_r, circle2_r, epsilon)) {
        // 圆重合
        return -1;
    }
    else if (fuzzy_equal(roundcast<T>(d), circle1_r + circle2_r, epsilon) || fuzzy_equal(jeabs(circle1_r - circle2_r), roundcast<T>(d), epsilon)) {
        // 相切
        jereal a = (jesqr(circle1_r) - jesqr(circle2_r) + jesqr(d)) / (2.0 * d);
        jereal ratio = jemul(a, d);
        inter_x1 = roundcast<T>(circle1_x + jemul(dx, ratio));
        inter_y1 = roundcast<T>(circle1_y + jemul(dy, ratio));
        return 1;
    }
    else {
        // 相交
        jereal a = (jesqr(circle1_r) - jesqr(circle2_r) + jesqr(d)) / (2.0 * d);
        jereal h = jesqrt(jesqr(circle1_r) - jesqr(a));
        jereal x2 = circle1_x + jemul(a, dx) / d;
        jereal y2 = circle1_y + jemul(a, dy) / d;
        inter_x1 = roundcast<T>(x2 + jemul(h, dy) / d);
        inter_y1 = roundcast<T>(y2 - jemul(h, dx) / d);
        inter_x2 = roundcast<T>(x2 - jemul(h, dy) / d);
        inter_y2 = roundcast<T>(y2 + jemul(h, dx) / d);
        return 2;
    }
}


template<typename T>
inline int intersection_arc2arc(const T& arc1_xc, const T& arc1_yc, const T& arc1_r, const T& arc1_start_rad, const T& arc1_end_rad,
                                 const T& arc2_xc, const T& arc2_yc, const T& arc2_r, const T& arc2_start_rad, const T& arc2_end_rad,
                                 T& inter_x1, T& inter_y1, T& inter_x2, T& inter_y2,
                                 const T& epsilon = Epsilon<T>()) noexcept
{
    T ix1, iy1, ix2, iy2;
    int ci_num = intersection_circle2circle<T>(arc1_xc, arc1_yc, arc1_r, arc2_xc, arc2_yc, arc2_r, ix1, iy1, ix2, iy2, epsilon);
    if (ci_num <= 0) return ci_num;
    int arc_inter_count = 0;
    if (ci_num >= 1) {
        T rad1 = vector_radian<T>(arc1_xc, arc1_yc, ix1, iy1);
        T rad2 = vector_radian<T>(arc2_xc, arc2_yc, ix1, iy1);
        if (is_radian_between(rad1, arc1_start_rad, arc1_end_rad, epsilon) && is_radian_between(rad2, arc2_start_rad, arc2_end_rad, epsilon)) {
            inter_x1 = ix1;
            inter_y1 = iy1;
            arc_inter_count++;
        }

    }
    if (ci_num >= 2) {
        T rad1 = vector_radian<T>(arc1_xc, arc1_yc, ix2, iy2);
        T rad2 = vector_radian<T>(arc2_xc, arc2_yc, ix2, iy2);
        if (is_radian_between(rad1, arc1_start_rad, arc1_end_rad, epsilon) && is_radian_between(rad2, arc2_start_rad, arc2_end_rad, epsilon)) {
            if (arc_inter_count > 0) {
                inter_x2 = ix2;
                inter_y2 = iy2;
            }
            else {
                inter_x1 = ix2;
                inter_y1 = iy2;
            }
            arc_inter_count++;
        }
    }
    return arc_inter_count;
}


/**
 * @brief 获取水平射线和直线段的交点
 */
template<typename T>
inline int intersection_rray2segment(const T &ray_x, const T &ray_y,
                                     const T &seg_sx, const T &seg_sy, const T &seg_ex, const T &seg_ey,
                                     T& inter_x1, T& inter_y1,
                                     const T& epsilon = Epsilon<T>())
{
    int seg_inter_count = 0;
    T dx = seg_sx - seg_ex;
    T dy = seg_sy - seg_ey;
    T smaxx = dx > 0 ? seg_sx : seg_ex;
    T sminy = dy < 0 ? seg_sy : seg_ey;
    T smaxy = dy > 0 ? seg_sy : seg_ey;

    if (ray_x > smaxx || ray_y < sminy || ray_y > smaxy) {
        return seg_inter_count;
    }

    if (dx == 0) {
        inter_x1 = smaxx;
        inter_y1 = ray_y;
        seg_inter_count++;
    } else {
        T seg_x = jediv(jemul(dx, (ray_y - seg_sy)), dy) + seg_sx;
        if (fuzzy_greaterequal(seg_x, ray_x, epsilon)) {
            inter_x1 = seg_x;
            inter_y1 = ray_y;
            seg_inter_count++;
        }
    }

    return seg_inter_count;
}


template<typename T>
inline bool intersection_rect2rect(const T &recta_x1, const T &recta_y1, const T &recta_x2, const T &recta_y2,
                                   const T &rectb_x1, const T &rectb_y1, const T &rectb_x2, const T &rectb_y2,
                                   T &inter_x1, T &inter_y1, T &inter_x2, T &inter_y2)
{
    T xmin = jemax(jemin(recta_x1, recta_x2), jemin(rectb_x1, rectb_x2));
    T xmax = jemin(jemax(recta_x1, recta_x2), jemax(rectb_x1, rectb_x2));
    T ymin = jemax(jemin(recta_y1, recta_y2), jemin(rectb_y1, rectb_y2));
    T ymax = jemin(jemax(recta_y1, recta_y2), jemax(rectb_y1, rectb_y2));

    if (xmin <= xmax && ymin <= ymax) {
        inter_x1 = xmin;
        inter_y1 = ymin;
        inter_x2 = xmax;
        inter_y2 = ymax;
        return true;
    }

    return false;
}


template<typename T>
inline bool intersects_segment2monarc(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                      const T& arc_xs, const T& arc_ys, const T& arc_xe, const T& arc_ye,
                                      const T& arc_xc, const T& arc_yc, const T& arc_r,
                                      bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (!intersects_rect2rect<T>(seg_xs, seg_ys, seg_xe, seg_ye, arc_xs, arc_ys, arc_xe, arc_ye, true, epsilon)) {
        return false;
    }

    T ix1, iy1, ix2, iy2;
    int inter_count = intersection_line2circle(seg_xs, seg_ys, seg_xe, seg_ye, arc_xc, arc_yc, arc_r, ix1, iy1, ix2, iy2, epsilon);
    if (inter_count <= 0) return false;


    T xmin = jemax(jemin(arc_xs, arc_xe), jemin(seg_xs, seg_xe));
    T xmax = jemin(jemax(arc_xs, arc_xe), jemax(seg_xs, seg_xe));
    T ymin = jemax(jemin(arc_ys, arc_ye), jemin(seg_ys, seg_ye));
    T ymax = jemin(jemax(arc_ys, arc_ye), jemax(seg_ys, seg_ye));


    if (inter_count == 1) {
        if (consider_touch == false) return false;
        if (contains_point_in_rect(ix1, iy1, xmin, ymin, xmax, ymax, true, epsilon)) {
            return true;
        }
    }
    else {
        if (contains_point_in_rect(ix1, iy1, xmin, ymin, xmax, ymax, consider_touch, epsilon)) {
            return true;
        }
        if (contains_point_in_rect(ix2, iy2, xmin, ymin, xmax, ymax, consider_touch, epsilon)) {
            return true;
        }
    }
}

template<typename T>
inline bool intersects_monarc2monarc(const T& arc1_xs, const T& arc1_ys, const T& arc1_xe, const T& arc1_ye,
                                     const T& arc1_xc, const T& arc1_yc, const T& arc1_r,
                                     const T& arc2_xs, const T& arc2_ys, const T& arc2_xe, const T& arc2_ye,
                                     const T& arc2_xc, const T& arc2_yc, const T& arc2_r,
                                     bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    T aminx1 = jemin(arc1_xs, arc1_xe);
    T aminy1 = jemin(arc1_ys, arc1_ye);
    T amaxx1 = jemax(arc1_xs, arc1_xe);
    T amaxy1 = jemax(arc1_ys, arc1_ye);
    T aminx2 = jemin(arc2_xs, arc2_xe);
    T aminy2 = jemin(arc2_ys, arc2_ye);
    T amaxx2 = jemax(arc2_xs, arc2_xe);
    T amaxy2 = jemax(arc2_ys, arc2_ye);
    if (jemax(aminx1, aminx2) > jemin(amaxx1, amaxx2)
        ||jemax(aminy1, aminy2) > jemin(amaxy1, amaxy2))
    {
        return false;
    }

    T minr = jemin(arc1_r, arc2_r);
    T maxr = jemax(arc1_r, arc2_r);
    T sumr = arc1_r + arc2_r;
    T distc = distance<T>(arc1_xc, arc1_yc, arc2_xc, arc2_yc);

    if (distc > 0 && distc <= sumr && distc + minr >= maxr) {
        // 相交, 求圆心连线和交点连线的交点 tc
        T dt1 = (jesqr(arc1_r) - jesqr(arc2_r) + jesqr(distc))/(2.0 * distc);
        T tuvx = jediv(arc2_xc - arc1_xc, distc);
        T tuvy = jediv(arc2_yc - arc1_yc, distc);
        T tvx = jemul(tuvx, dt1);
        T tvy = jemul(tuvy, dt1);
        T tcx = arc1_xc + tvx;
        T tcy = arc1_yc + tvy;
        if (fuzzy_equal(dt1, arc1_r, epsilon)) {
            // 相切
            if (!consider_touch)
                return false;
            if (consider_touch
                && fuzzy_greaterequal(tcx, aminx1, epsilon) && fuzzy_lessequal(tcx, amaxx1, epsilon)
                && fuzzy_greaterequal(tcy, aminy1, epsilon) && fuzzy_lessequal(tcy, amaxy1, epsilon)
                && fuzzy_greaterequal(tcx, aminx2, epsilon) && fuzzy_lessequal(tcx, amaxx2, epsilon)
                && fuzzy_greaterequal(tcy, aminy2, epsilon) && fuzzy_lessequal(tcy, amaxy2, epsilon))
            {
                return true;
            }
            else if (!consider_touch
                       && tcx > aminx1 && tcx < amaxx1
                       && tcy > aminy1 && tcy < amaxy1
                       && tcx > aminx2 && tcx < amaxx2
                       && tcy > aminy2 && tcy < amaxy2)
            {
                return true;
            }
        }
        else {
            // 两个交点, 圆上交点到 tc 的距离
            T dpttc = jesqrt<T>(jesqr(arc1_r) - jesqr(dt1));
            // 求出交点
            T px1 = roundcast<T>(tcx + jemul(-tuvy, dpttc));
            T py1 = roundcast<T>(tcy + jemul(tuvx, dpttc));
            if (consider_touch
                && fuzzy_greaterequal(px1, aminx1, epsilon) && fuzzy_lessequal(px1, amaxx1, epsilon)
                && fuzzy_greaterequal(py1, aminy1, epsilon) && fuzzy_lessequal(py1, amaxy1, epsilon)
                && fuzzy_greaterequal(px1, aminx2, epsilon) && fuzzy_lessequal(px1, amaxx2, epsilon)
                && fuzzy_greaterequal(py1, aminy2, epsilon) && fuzzy_lessequal(py1, amaxy2, epsilon))
            {
                return true;
            }
            else if (!consider_touch
                       && px1 > aminx1 && px1 < amaxx1
                       && py1 > aminy1 && py1 < amaxy1
                       && px1 > aminx2 && px1 < amaxx2
                       && py1 > aminy2 && py1 < amaxy2)
            {
                return true;
            }

            T px2 = roundcast<T>(tcx + jemul(tuvy, dpttc));
            T py2 = roundcast<T>(tcy + jemul(-tuvx, dpttc));
            if (consider_touch
                && fuzzy_greaterequal(px2, aminx1, epsilon) && fuzzy_lessequal(px2, amaxx1, epsilon)
                && fuzzy_greaterequal(py2, aminy1, epsilon) && fuzzy_lessequal(py2, amaxy1, epsilon)
                && fuzzy_greaterequal(px2, aminx2, epsilon) && fuzzy_lessequal(px2, amaxx2, epsilon)
                && fuzzy_greaterequal(py2, aminy2, epsilon) && fuzzy_lessequal(py2, amaxy2, epsilon))
            {
                return true;
            }
            else if (!consider_touch
                       && px2 > aminx1 && px2 < amaxx1
                       && py2 > aminy1 && py2 < amaxy1
                       && px2 > aminx2 && px2 < amaxx2
                       && py2 > aminy2 && py2 < amaxy2)
            {
                return true;
            }
        }
    }
    return false;
}


template<typename T>
inline bool intersects_hsegment2monrect(const T& hseg_y, const T& hseg_xmin, const T& hseg_xmax,
                                        const T& rect_xmin, const T& rect_ymin, const T& rect_xmax, const T& rect_ymax,
                                        bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (!is_value_between(hseg_y, rect_ymin, rect_ymax, consider_touch, epsilon)) {
        return false;
    }
    return is_interval_cross(hseg_xmin, hseg_xmax, rect_xmin, rect_xmax, consider_touch, epsilon);
}

template<typename T>
inline bool intersects_vsegment2monrect(const T& vseg_x, const T& vseg_ymin, const T& vseg_ymax,
                                        const T& rect_xmin, const T& rect_ymin, const T& rect_xmax, const T& rect_ymax,
                                        bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (!is_value_between(vseg_x, rect_xmin, rect_xmax, consider_touch, epsilon)) {
        return false;
    }
    return is_interval_cross(vseg_ymin, vseg_ymax, rect_ymin, rect_ymax, consider_touch, epsilon);
}

template<typename T>
inline bool intersects_rect2monarc(const T& rect_x1, const T& rect_y1, const T& rect_x2, const T& rect_y2,
                                   const T& arc_xs, const T& arc_ys, const T& arc_xe, const T& arc_ye,
                                   const T& arc_xc, const T& arc_yc, const T& arc_r,
                                   bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    T axmin = jemin(arc_xs, arc_xe);
    T axmax = jemax(arc_xs, arc_xe);
    T aymin = jemin(arc_ys, arc_ye);
    T aymax = jemax(arc_ys, arc_ye);
    T rxmin = jemin(rect_x1, rect_x2);
    T rxmax = jemax(rect_x1, rect_x2);
    T rymin = jemin(rect_y1, rect_y2);
    T rymax = jemax(rect_y1, rect_y2);
    if (contains_monrect_in_monrect(axmin, aymin, axmax, aymax, rxmin, rymin, rxmax, rymax, consider_touch, epsilon)) {
        return true;
    }

    if (contains_point_in_monrect(arc_xs, arc_ys, rxmin, rymin, rxmax, rymax, consider_touch, epsilon)) {
        return true;
    }
    if (contains_point_in_monrect(arc_xe, arc_ye, rxmin, rymin, rxmax, rymax, consider_touch, epsilon)) {
        return true;
    }

    T inter_minx = jemax(axmin, rxmin);
    T inter_maxx = jemin(axmax, rxmax);
    T inter_miny = jemax(aymin, rymin);
    T inter_maxy = jemin(aymax, rymax);
    if (intersects_hsegment2monrect(rymax, rxmin, rxmax, axmin, aymin, axmax, aymax, consider_touch, epsilon)) {
        T dx = jesqrt<T>(jesqr(arc_r) - jesqr(rymax - arc_yc));
        T x1 = arc_xc + dx;
        T x2 = arc_xc - dx;
        return (is_value_between(x1, inter_minx, inter_maxx, consider_touch, epsilon)
                || is_value_between(x2, inter_minx, inter_maxx, consider_touch, epsilon));
    }
    if (intersects_hsegment2monrect(rymin, rxmin, rxmax, axmin, aymin, axmax, aymax, consider_touch, epsilon)) {
        T dx = jesqrt<T>(jesqr(arc_r) - jesqr(rymin - arc_yc));
        T x1 = arc_xc + dx;
        T x2 = arc_xc - dx;
        return (is_value_between(x1, inter_minx, inter_maxx, consider_touch, epsilon)
                || is_value_between(x2, inter_minx, inter_maxx, consider_touch, epsilon));
    }
    if (intersects_vsegment2monrect(rxmax, rymin, rymax, axmin, aymin, axmax, aymax, consider_touch, epsilon)) {
        T dy = jesqrt<T>(jesqr(arc_r) - jesqr(rxmax - arc_xc));
        T y1 = arc_yc + dy;
        T y2 = arc_yc - dy;
        return (is_value_between(y1, inter_miny, inter_maxy, consider_touch, epsilon)
                || is_value_between(y2, inter_miny, inter_maxy, consider_touch, epsilon));
    }
    if (intersects_vsegment2monrect(rxmin, rymin, rymax, axmin, aymin, axmax, aymax, consider_touch, epsilon)) {
        T dy = jesqrt<T>(jesqr(arc_r) - jesqr(rxmin - arc_xc));
        T y1 = arc_yc + dy;
        T y2 = arc_yc - dy;
        return (is_value_between(y1, inter_miny, inter_maxy, consider_touch, epsilon)
                || is_value_between(y2, inter_miny, inter_maxy, consider_touch, epsilon));
    }
    return false;
}



template<typename T>
inline bool intersects_segment2arc(const T& seg_xs, const T& seg_ys, const T& seg_xe, const T& seg_ye,
                                   const T& arc_xs, const T& arc_ys, const T& arc_xe, const T& arc_ye,
                                   const T& arc_xc, const T& arc_yc, const T& arc_r,
                                   const jereal &arc_start_rad, const jereal &arc_end_rad,
                                   bool consider_touch = true, const T& epsilon = Epsilon<T>()) noexcept
{
    if (!intersects_rect2rect<T>(seg_xs, seg_ys, seg_xe, seg_ye, arc_xs, arc_ys, arc_xe, arc_ye, true, epsilon)) {
        return false;
    }

    T ix1, iy1, ix2, iy2;
    int inter_count = intersection_line2circle(seg_xs, seg_ys, seg_xe, seg_ye, arc_xc, arc_yc, arc_r, ix1, iy1, ix2, iy2, epsilon);
    if (inter_count <= 0) return false;
    if (inter_count == 1) {
        if (consider_touch == false) return false;
        jereal rad = vector_radian<T>(arc_xc, arc_yc, ix1, iy1);
        if (is_radian_between(rad, arc_start_rad, arc_end_rad, epsilon)
            && is_point_on_segment(ix1, iy1, seg_xs, seg_ys, seg_xe, seg_ye, epsilon))
        {
            return true;
        }
    }
    else {
        jereal rad1 = vector_radian<T>(arc_xc, arc_yc, ix1, iy1);
        if (is_radian_between(rad1, arc_start_rad, arc_end_rad, epsilon) && is_point_on_segment(ix1, iy1, seg_xs, seg_ys, seg_xe, seg_ye, epsilon)) {
            if (consider_touch == true) {
                return true;
            } else {
                if (is_point_equal(ix1, iy1, seg_xs, seg_ys) || is_point_equal(ix1, iy1, seg_xe, seg_ye)
                    || is_point_equal(ix1, iy1, arc_xs, arc_ys) || is_point_equal(ix1, iy1, arc_xe, arc_ye)) {
                    return false;
                } else {
                    return true;
                }
            }
        }
        jereal rad2 = vector_radian<T>(arc_xc, arc_yc, ix2, iy2);
        if (is_radian_between(rad2, arc_start_rad, arc_end_rad, epsilon) && is_point_on_segment(ix2, iy2, seg_xs, seg_ys, seg_xe, seg_ye, epsilon)) {
            if (consider_touch == true) {
                return true;
            } else {
                if (is_point_equal(ix2, iy2, seg_xs, seg_ys) || is_point_equal(ix2, iy2, seg_xe, seg_ye)
                    || is_point_equal(ix2, iy2, arc_xs, arc_ys) || is_point_equal(ix2, iy2, arc_xe, arc_ye)) {
                    return false;
                } else {
                    return true;
                }
            }
        }
    }
}




/**
 * @brief 点到线段的最近点
 * @param px,py : 点坐标
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段起点和终点
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T>
inline void closest_point2segment(const T &px, const T &py,
                                  const T &seg_xs, const T &seg_ys, const T &seg_xe, const T &seg_ye,
                                  T &closest_x, T &closest_y) noexcept
{
    const T vx = seg_xe - seg_xs;
    const T vy = seg_ye - seg_ys;
    const T wx = px - seg_xs;
    const T wy = py - seg_ys;

    const jereal c1 = dot_product<vx, vy, wx, wy>;
    if (c1 <= 0) {
        closest_x = seg_xs;
        closest_y = seg_ys;
        return;
    }

    const jereal c2 = dot_product<vx, vy, vx, vy>;
    if (c2 <= c1) {
        closest_x = seg_xe;
        closest_y = seg_ye;
        return;
    }
    const jereal t = c1 / c2;
    closest_x = roundcast<T, jereal>(seg_xs + t * vx);
    closest_y = roundcast<T, jereal>(seg_ys + t * vy);
}

/**
 * @brief 点到射线的最近点
 * @param px, py : 点坐标
 * @param ray_xo, ray_yo : 射线起点
 * @param ray_xd, ray_yd : 射线方向
 * @param closest_x, closest_y : 返回的最近点坐标
 * @return 如果点在射线起点前,则没有最近点, 返回0, 否则返回1
 */
template<typename T>
inline void closest_point2ray(const T &px, const T &py,
                              const T &ray_xo, const T &ray_yo, const T &ray_xd,
                              const T &ray_yd, T &closest_x, T &closest_y) noexcept
{
    const jereal t = dot_product(ray_xd, ray_yd, px - ray_xo, py - ray_yo);
    if (t <= 0) {
        closest_x = ray_xo;
        closest_y = ray_yo;
        return;
    }
    else {
        closest_x = roundcast<T, jereal>(ray_xo + t * ray_xd);
        closest_y = roundcast<T, jereal>(ray_yo + t * ray_yd);
        return;
    }
}

/**
 * @brief 点到矩形的最近点
 * @param px, py : 点坐标
 * @param rect_x1, rect_y1, rect_x2, rect_y2 : 矩形的对角线两个点
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T>
inline void closest_point2rect(const T &px, const T &py,
                               const T &rect_x1, const T &rect_y1, const T &rect_x2, const T &rect_y2,
                               T &closest_x, T &closest_y) noexcept
{
    if (px < jemin(rect_x1, rect_x2))
        closest_x = jemin(rect_x1, rect_x2);
    else if (px > jemax(rect_x1, rect_x2))
        closest_x = jemax(rect_x1, rect_x2);
    else
        closest_x = px;

    if (py < jemin(rect_y1, rect_y2))
        closest_y = jemin(rect_y1, rect_y2);
    else if (py > jemax(rect_y1, rect_y2))
        closest_y = jemax(rect_y1, rect_y2);
    else
        closest_y = py;
}

/**
 * @brief 点到圆的最近点, 如果在圆内则返回点坐标
 * @param px, py : 点坐标
 * @param circle_x, circle_y : 圆心坐标
 * @param circle_r : 圆半径
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T >
inline void closest_point2circle(const T &px, const T &py,
                                 const T &circle_x, const T &circle_y, const T &circle_r,
                                 T &closest_x, T &closest_y) noexcept
{
    const T dx = px - circle_x;
    const T dy = py - circle_y;
    if ((jesqr(dx) + jesqr(dy)) <= jesqr(circle_r)) {
        closest_x = px;
        closest_y = py;
    }
    else {
        const jereal ratio = jediv(circle_r, jehypot(dx, dy));
        closest_x = roundcast<T, jereal>(circle_x + dx * ratio);
        closest_y = roundcast<T, jereal>(circle_y + dy * ratio);
    }
}

/**
 * @brief 线段到圆的最近点
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段起点和终点
 * @param circle_x, circle_y, circle_r : 圆心和半径
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T>
inline void closest_segment2circle(const T &seg_xs, const T &seg_ys, const T &seg_xe, const T &seg_ye,
                                   const T &circle_x, const T &circle_y, const T &circle_r,
                                   T &closest_x, T &closest_y) noexcept
{
    T nx = 0.0;
    T ny = 0.0;
    closest_point2segment<T>(circle_x, circle_y, seg_xs, seg_ys, seg_xe, seg_ye, nx, ny);
    const jereal ratio = jediv(circle_r, jehypot(nx - circle_x, ny - circle_y));
    closest_x = roundcast<T, jereal>(circle_x + static_cast<jereal>((nx - circle_x) * ratio));
    closest_y = roundcast<T, jereal>(circle_y + static_cast<jereal>((ny - circle_y) * ratio));
}

/**
 * @brief 直线到圆的最近点
 * @param line_x1, line_y1, line_x2, line_y2 : 直线起点和终点
 * @param circle_x, circle_y, circle_r : 圆心和半径
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T>
inline void closeset_line2circle(const T &line_x1, const T &line_y1, const T &line_x2, const T &line_y2,
                                 const T &circle_x, const T &circle_y, const T &circle_r,
                                 T &closest_x, T &closest_y) noexcept
{
    T nx = 0.0;
    T ny = 0.0;
    closest_point2line<T>(circle_x, circle_y, line_x1, line_y1, line_x2, line_y2, nx, ny);
    const jereal ratio = jediv(circle_r, jehypot(nx - circle_x, ny - circle_y));
    closest_x = roundcast<T, jereal>(circle_x + static_cast<jereal>((nx - circle_x) * ratio));
    closest_y = roundcast<T, jereal>(circle_y + static_cast<jereal>((ny - circle_y) * ratio));
}



/**
 * @brief 圆到圆的最近点
 * @param circle1_x, circle1_y, circle1_r : 圆1的圆心和半径
 * @param circle2_x, circle2_y, circle2_r : 圆2的圆心和半径
 * @param closest1_x, closest1_y : 返回的最近点1坐标
 * @param closest2_x, closest2_y : 返回的最近点2坐标
 * @return 返回最近点的个数,(0:表示两圆相交, 1:表示两圆相切, 2:表示两圆不相交)
 */
template<typename T >
inline int closest_circle2circle(const T &circle1_x, const T &circle1_y, const T &circle1_r,
                                 const T &circle2_x, const T &circle2_y, const T &circle2_r,
                                 T &closest1_x, T &closest1_y,T &closest2_x, T& closest2_y) noexcept
{
    const T dx = circle2_x - circle1_x;
    const T dy = circle2_y - circle1_y;
    const T d = jehypot(dx, dy);
    if (fuzzy_equal(d, circle1_r + circle2_r)) {
        // 两圆相切, 则切点为最近点
        closest1_x = roundcast<T, jereal>(circle1_x + dx * static_cast<jereal>(circle1_r) / d);
        closest1_y = roundcast<T, jereal>(circle1_y + dy * static_cast<jereal>(circle1_r) / d);
        return 1;
    }
    else if (d > circle1_r + circle2_r) {
        // 两圆相离, 则中心连线分别与两圆的交点为最近点
        const jereal ratio = jediv(circle1_r, d);
        closest1_x = roundcast<T, jereal>(circle1_x + dx * ratio);
        closest1_y = roundcast<T, jereal>(circle1_y + dy * ratio);

        ratio = jediv(circle2_r, d);
        closest2_x = roundcast<T, jereal>(circle2_x - dx * ratio);
        closest2_y = roundcast<T, jereal>(circle2_y - dy * ratio);
        return 2;
    }
    else {
        // 两圆相交, 则没有最近点
        return 0;
    }
}


/**
 * @brief 计算点到圆弧的最近点
 * @param px, py : 点坐标
 * @param arc_xs, arc_ys, arc_xe, arc_ye : 圆弧的起点和终点
 * @param arc_xc, arc_yc, arc_r : 圆弧的圆心,半径
 * @param arc_sa_rad, arc_ea_rad : 圆弧的起始角度和结束角度
 * @param closest_x, closest_y : 返回的最近点坐标
 */
template<typename T>
inline void closest_point2arc(const T &px, const T &py,
                              const T &arc_xs, const T &arc_ys, const T &arc_xe, const T &arc_ye,
                              const T &arc_xc, const T &arc_yc, const T &arc_r,
                              const jereal &arc_sa_rad, const jereal &arc_ea_rad,
                              T &closest_x, T &closest_y) noexcept
{
    const T dx = px - arc_xc;
    const T dy = py - arc_yc;
    // 如果点在圆心,则最近点为圆弧的起点
    if (fuzzy_iszero(dx) && fuzzy_iszero(dy)) {
        closest_x = arc_xs;
        closest_y = arc_ys;
        return ;
    }

    const jereal rad = std::atan2(dx, dy);
    if (is_radian_between(rad, arc_sa_rad, arc_ea_rad)) {
        const jereal ratio = jediv(arc_r, jehypot(dx, dy));
        closest_x = roundcast<T, jereal>(arc_xc + dx * ratio);
        closest_y = roundcast<T, jereal>(arc_yc + dy * ratio);
        return;
    } else {
        auto ds = sqrdist<T>(px, py, arc_xs, arc_ys);
        auto de = sqrdist<T>(px, py, arc_xe, arc_ye);
        if (ds < de) {
            closest_x = arc_xs;
            closest_y = arc_ys;
        } else {
            closest_x = arc_xe;
            closest_y = arc_ye;
        }
    }
}

/**
 * @brief 点到点的平方距离
 * @param p1xs,p1ys : 点1坐标
 * @param p2xs,p2ys : 点2坐标
 * @return
 */
template<typename T>
inline jereal sqrdist_point2point(const T &p1xs, const T &p1ys,
                                  const T &p2xs, const T& p2ys) noexcept
{
    return sqrdist<T>(p1xs, p1ys, p2xs, p2ys);
}

/**
 * @brief 点到点的距离
 * @param p1xs,p1ys : 点1坐标
 * @param p2xs,p2ys : 点2坐标
 * @return
 */
template<typename T>
inline T distance_point2point(const T &p1xs, const T &p1ys,
                              const T &p2xs, const T& p2ys) noexcept
{
    return distance<T>(p1xs, p1ys, p2xs, p2ys);
}



/**
 * @brief 点到直线的平方距离
 * @param px, py : 点坐标
 * @param line_x1, line_y1, line_x2, line_y2 : 直线起点和终点
 * @return
 */
template<typename T>
inline jereal sqrdist_point2line(const T &px, const T &py,
                                 const T &line_x1, const T &line_y1, const T &line_x2, const T &line_y2) noexcept
{
    T closest_x, closest_y;
    closest_point2line<T>(px, py, line_x1, line_y1, line_x2, line_y2, closest_x, closest_y);
    return sqrdist_point2point<T>(px, py, closest_x, closest_y);
}

/**
 * @brief 点到直线的距离
 * @param px, py : 点坐标
 * @param line_x1, line_y1, line_x2, line_y2 : 直线起点和终点
 * @return
 */
template<typename T>
inline T distance_point2line(const T &px, const T &py,
                             const T &line_x1, const T &line_y1, const T &line_x2, const T &line_y2) noexcept
{
    T closest_x, closest_y;
    closest_point2line<T>(px, py, line_x1, line_y1, line_x2, line_y2, closest_x, closest_y);
    return distance_point2point<T>(px, py, closest_x, closest_y);
}

/**
 * @brief 点到线段的平方距离
 * @param px, py : 点坐标
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段起点和终点
 * @return
 */
template<typename T>
inline jereal sqrdist_point2segment(const T &px, const T &py,
                                    const T &seg_xs, const T &seg_ys, const T &seg_xe, const T &seg_ye) noexcept
{
    T closest_x, closest_y;
    closest_point2segment<T>(px, py, seg_xs, seg_ys, seg_xe, seg_ye, closest_x, closest_y);
    return sqrdist_point2point<T>(px, py, closest_x, closest_y);
}

/**
 * @brief 点到线段的距离
 * @param px, py : 点坐标
 * @param seg_xs, seg_ys, seg_xe, seg_ye : 线段起点和终点
 * @return
 */
template<typename T>
inline T distance_point2segment(const T &px, const T &py,
                                const T &seg_xs, const T &seg_ys, const T &seg_xe, const T &seg_ye) noexcept
{
    T closest_x, closest_y;
    closest_point2segment<T>(px, py, seg_xs, seg_ys, seg_xe, seg_ye, closest_x, closest_y);
    return distance_point2point<T>(px, py, closest_x, closest_y);
}

/**
 * @brief 点到射线的平方距离
 * @param px, py : 点坐标
 * @param ray_xo, ray_yo : 射线起点
 * @param ray_xd, ray_yd : 射线方向
 * @return
 */
template<typename T>
inline jereal sqrdist_point2ray(const T &px, const T &py,
                                const T &ray_xo, const T &ray_yo, const T &ray_xd, const T &ray_yd) noexcept
{
    T closest_x, closest_y;
    closest_point2ray<T>(px, py, ray_xo, ray_yo, ray_xd, ray_yd, closest_x, closest_y);
    return sqrdist_point2point<T>(px, py, closest_x, closest_y);
}

/**
 * @brief 点到射线的距离
 * @param px, py : 点坐标
 * @param ray_xo, ray_yo : 射线起点
 * @param ray_xd, ray_yd : 射线方向
 * @return
 */
template<typename T>
inline T distance_point2ray(const T &px, const T &py,
                            const T &ray_xo, const T &ray_yo, const T &ray_xd, const T &ray_yd) noexcept
{
    T closest_x, closest_y;
    closest_point2ray<T>(px, py, ray_xo, ray_yo, ray_xd, ray_yd, closest_x, closest_y);
    return distance_point2point<T>(px, py, closest_x, closest_y);
}


/**
 * @brief 点到矩形的平方距离
 * @param px, py : 点坐标
 * @param rect_x1, rect_y1, rect_x2, rect_y2 : 矩形的对角线两个点
 * @return
 */
template<typename T>
inline jereal sqrdist_point2rect(const T &px, const T &py,
                                 const T &rect_x1, const T &rect_y1, const T &rect_x2, const T &rect_y2) noexcept
{
    T closest_x, closest_y;
    closest_point2rect<T>(px, py, rect_x1, rect_y1, rect_x2, rect_y2, closest_x, closest_y);
    return sqrdist_point2point<T>(px, py, closest_x, closest_y);
}

/**
 * @brief 点到矩形的距离
 * @param px, py : 点坐标
 * @param rect_x1, rect_y1, rect_x2, rect_y2 : 矩形的对角线两个点
 * @return
 */
template<typename T>
inline T distance_point2rect(const T &px, const T &py,
                             const T &rect_x1, const T &rect_y1, const T &rect_x2, const T &rect_y2) noexcept
{
    T closest_x, closest_y;
    closest_point2rect<T>(px, py, rect_x1, rect_y1, rect_x2, rect_y2, closest_x, closest_y);
    return distance_point2point<T>(px, py, closest_x, closest_y);
}


/**
 * @brief 直线到直线的平方距离
 * @param line1_xs, line1_ys, line1_xe, line1_ye : 直线1起点和终点
 * @param line2_xs, line2_ys, line2_xe, line2_ye : 直线2起点和终点
 */
template<typename T=jereal>
inline jereal sqrdist_line2line(const T &line1_xs, const T &line1_ys, const T &line1_xe, const T &line1_ye,
                                const T &line2_xs, const T &line2_ys, const T &line2_xe, const T &line2_ye) noexcept
{
    const jereal ux = line1_xe - line1_xs;
    const jereal uy = line1_ye - line1_ys;
    const jereal vx = line2_xe - line2_xs;
    const jereal vy = line2_ye - line2_ys;

    if (fuzzy_notequal<T>(ux * vy,  uy * vx)) {
        return 0.0;
    }

    const jereal wx = line1_xs - line2_xs;
    const jereal wy = line1_ys - line2_ys;

    const jereal a  = (ux * ux + uy * uy);
    const jereal b  = (ux * vx + uy * vy);
    const jereal c  = (vx * vx + vy * vy);
    const jereal d  = (ux * wx + uy * wy);
    const jereal e  = (vx * wx + vy * wy);
    const jereal dt = a * c - b * b;

    jereal sc = 0.0;
    jereal tc = 0.0;

    if (fuzzy_equal(dt, T(0.0))) {
        sc = T(0.0);

        if (b > c)
            tc = d / b;
        else
            tc = e / c;
    }
    else {
        sc = (b * e - c * d) / dt;
        tc = (a * e - b * d) / dt;
    }

    const jereal dx = wx + (sc * ux) - (tc * vx);
    const jereal dy = wy + (sc * uy) - (tc * vy);

    return dx * dx + dy * dy;
}

/**
 * @brief 直线到直线的距离
 * @param line1_xs, line1_ys, line1_xe, line1_ye : 直线1起点和终点
 * @param line2_xs, line2_ys, line2_xe, line2_ye : 直线2起点和终点
 */
template<typename T>
inline T distance_line2line(const T &line1_xs, const T &line1_ys, const T &line1_xe, const T &line1_ye,
                            const T &line2_xs, const T &line2_ys, const T &line2_xe, const T &line2_ye) noexcept
{
    return jesqrt<T>(sqrdist_line2line<T>(line1_xs, line1_ys, line1_xe, line1_ye, line2_xs, line2_ys, line2_xe, line2_ye));
}

/**
 * @brief 线段到线段的平方距离
 * @param seg1_xs, seg1_ys, seg1_xe, seg1_ye : 线段1起点和终点
 * @param seg2_xs, seg2_ys, seg2_xe, seg2_ye : 线段2起点和终点
 */
template<typename T>
inline jereal sqrdist_segment2segment(const T &seg1_xs, const T &seg1_ys, const T &seg1_xe, const T &seg1_ye,
                                      const T &seg2_xs, const T &seg2_ys, const T &seg2_xe, const T &seg2_ye) noexcept
{
    const jereal ux = seg1_xe - seg1_xs;
    const jereal uy = seg1_ye - seg1_ys;

    const jereal vx = seg2_xe - seg2_xs;
    const jereal vy = seg2_ye - seg2_ys;

    const jereal wx = seg1_xs - seg2_xs;
    const jereal wy = seg1_ys - seg2_ys;

    const jereal a = (ux * ux + uy * uy);
    const jereal b = (ux * vx + uy * vy);
    const jereal c = (vx * vx + vy * vy);
    const jereal d = (ux * wx + uy * wy);
    const jereal e = (vx * wx + vy * wy);
    const jereal dt = a * c - b * b;

    jereal sc = 0.0;
    jereal sn = 0.0;
    jereal tc = 0.0;
    jereal tn = 0.0;
    jereal sd = dt;
    jereal td = dt;

    if (fuzzy_equal(dt, 0.0)) {
        sn = 0.0;
        sd = 1.0;
        tn = e;
        td = c;
    } else {
        sn = (b * e - c * d);
        tn = (a * e - b * d);

        if (sn < 0.0) {
            sn = 0.0;
            tn = e;
            td = c;
        } else if (sn > sd) {
            sn = sd;
            tn = e + b;
            td = c;
        }
    }

    if (tn < 0.0) {
        tn = 0.0;

        if (-d < 0.0)
            sn = 0.0;
        else if (-d > a)
            sn = sd;
        else {
            sn = -d;
            sd = a;
        }
    } else if (tn > td) {
        tn = td;

        if ((-d + b) < 0.0)
            sn = 0.0;
        else if ((-d + b) > a)
            sn = sd;
        else {
            sn = (-d + b);
            sd = a;
        }
    }

    if (fuzzy_equal(sn, 0.0))
        sc = 0.0;
    else
        sc = sn / sd;

    if (fuzzy_equal(tn, 0.0))
        tc = 0.0;
    else
        tc = tn / td;

    const jereal dx = wx + (sc * ux) - (tc * vx);
    const jereal dy = wy + (sc * uy) - (tc * vy);

    return dx * dx + dy * dy;
}


/**
 * @brief 线段到线段的距离
 * @param seg1_xs, seg1_ys, seg1_xe, seg1_ye : 线段1起点和终点
 * @param seg2_xs, seg2_ys, seg2_xe, seg2_ye : 线段2起点和终点
 */
template<typename T>
inline T distance_segment2segment(const T &seg1_xs, const T &seg1_ys, const T &seg1_xe, const T &seg1_ye,
                                  const T &seg2_xs, const T &seg2_ys, const T &seg2_xe, const T &seg2_ye) noexcept
{
    return jesqrt<T>(sqrdist_segment2segment<T>(seg1_xs, seg1_ys, seg1_xe, seg1_ye, seg2_xs, seg2_ys, seg2_xe, seg2_ye));
}

template<typename T>
std::vector<std::tuple<T,T>> arc_monotone_points(const T& xs, const T& ys, const T& xe, const T& ye, const T& xc, const T& yc, bool cw ) {
    T radius = jehypot<T>(xs - xc, ys - yc);
    jereal sa = std::atan2(ys - yc, xs - xc) * D180_DIV_PI;
    jereal ea = std::atan2(ye - yc, xe - xc) * D180_DIV_PI;
    sa = normalize_degree(sa);
    ea = normalize_degree(ea);
    std::vector<jereal> angles;
    if (sa >= 0 && sa < 90)
        angles = {90, 180, 270, 0};
    else if (sa >= 90 && sa < 180)
        angles = {180, 270, 0, 90};
    else if (sa >= 180 && sa < 270)
        angles = {270, 0, 90, 180};
    else if (sa >= 270 && sa < 360)
        angles = {0, 90, 180, 270};
    if (cw)
        std::reverse(angles.begin(), angles.end());
    std::vector<std::tuple<T, T>> points;
    for (auto angle : angles) {
        bool ok = false;
        if (cw) {
            ok = (sa >= ea) ? (angle < sa && angle > ea) : (angle < sa || angle > ea);
        } else {
            ok = (sa <= ea) ? (angle > sa && angle < ea) : (angle > sa || angle < ea);
        }
        if (ok) {
            if (angle == 0) {
                points.push_back(std::tuple<T,T>(xc + radius, yc));
            } else if (angle == 90) {
                points.push_back(std::tuple<T,T>(xc, yc + radius));
            } else if (angle == 180) {
                points.push_back(std::tuple<T,T>(xc - radius, yc));
            } else if (angle == 270) {
                points.push_back(std::tuple<T,T>(xc, yc - radius));
            }
        }
    }
    return points;
}


template<typename T>
std::vector<std::tuple<T,T>> arc_bound_points(const T& xs, const T& ys, const T& xe, const T& ye, const T& xc, const T& yc, bool cw ) {
    auto points = arc_monotone_points(xs, ys, xe, ye, xc, yc, cw);
    points.push_back(std::tuple<T,T>(xs, ys));
    points.push_back(std::tuple<T,T>(xe, ye));
    return points;
}


template<typename T>
std::vector<std::tuple<T, T>> arc_bound_points(const T &xs, const T &ys, const T &xe, const T &ye, jereal angle)
{
    bool cw = angle < 0;
    T xc, yc;
    arc_center_from_angle(xs, ys, xe, ye, angle, xc, yc);
    return arc_bound_points(xs, ys, xe, ye, xc, yc, cw);
}

} // namespace


#endif // JEGEO_MATH_HPP
